Introducción al análsis de imagenes satelitales

Este video exploraremos cómo importar y graficar un ráster multibanda en R. Usaremos un subset de una escena de Landsat 9 producto de surferce reflectance, perteneciente al bosque Chaqueño.

Primero aprenderemos a leer un raster, y visualizarrlo. Finalmente, calcularemos distintos indices de vegetación.

A continuación, te mostramos una comparativa de las bandas de ambos satélites (Landsat 9):

Bandas

Librerias a utilizar

# load the raster, sp, and rgdal packages
rm(list=ls())
library(raster)
## Loading required package: sp
library(sp)
library(rgdal)
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.2, released 2022/03/08
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(ggplot2)
# load raster in an R object
setwd("/Users/veronica/Documents/Maestria/modelado_espacial/modelado_espacial/trabajo_final/data")
# Blue
imagen <- raster('L9_chaco.tif')
imagen_stack = stack("./L9_chaco.tif")
summary(imagen)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (9.86% of all cells)
##           L9_chaco
## Min.    -0.0571925
## 1st Qu.  0.0179375
## Median   0.0197800
## 3rd Qu.  0.0228050
## Max.     0.4840350
## NA's     0.0000000
# create a grayscale color palette to use for the image.
grayscale_colors <- gray.colors(100,            # number of different color levels 
                                start = 0.0,    # how black (0) to go
                                end = 1.0,      # how white (1) to go
                                gamma = 2.2,    # correction between how a digital 
                                # camera sees the world and how human eyes see it
                                alpha = NULL)   #Null=colors are not transparent

# Plot band 1

# plot all three bands separately
plot(imagen_stack, 
     col=grayscale_colors)

# view attributes: Check out dimension, CRS, resolution, values attributes, and 
# band.
imagen_stack
## class      : RasterStack 
## dimensions : 789, 1285, 1013865, 8  (nrow, ncol, ncell, nlayers)
## resolution : 0.0002694946, 0.0002694946  (x, y)
## extent     : -62.29071, -61.94441, -25.94532, -25.73269  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      : SR_B1, SR_B2, SR_B3, SR_B4, SR_B5, SR_B6, SR_B7, ST_B10
# Can specify which band we want to read in
banda_2 <- 
  raster(paste0("L9_chaco.tif"), 
           band = 2)

# plot band 2
plot(banda_2,
     col=grayscale_colors, # we already created this palette & can use it again
     axes=FALSE,
     main="RGB Imagery- SR_B2-Green")

Vamos a analizar e interpretar el contenido de información de las bandas de las imágenes ópticas multiespectrales. Familiarizarnos con las diferentes opciones de combinación de bandas y la obtención de imágenes de color.

# view histogram of all 3 bands
hist(imagen_stack,
     maxpixels=ncell(imagen_stack))

Observen que la imagen de arriba se ve oscura. A continuación, utilizaremos el argumento de estiramiento (strech) integrado en la función plot_rgb (). El argumento str_clip le permite especificar la cantidad de partes de la distribución de datos se desea recortar. Cuanto mayor sea el número, más se estirarán o se aclararán los datos.

#landsatRGB <- stack(b4, b3, b2)
plotRGB(imagen_stack, 
        r = 'SR_B4', g = 'SR_B3', b = 'SR_B2',axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")

Las tonalidades más habituales en una composición en falso color son:

#landsatFCC <- stack(b5, b4, b3)
plotRGB(imagen_stack, 
        r = 'SR_B5', g = 'SR_B4', b = 'SR_B3', axes=TRUE, stretch="lin", main="Landsat False Color Composite")

## Indices de vegetación

Los índices son mediciones empíricas de las propiedades de la superficie terrestre; constituyen valores adimensionales basados en medidas de radiación y calculados a partir de la combinación de bandas espectrales.

El objetivo principal de los índices es formular una medida sintética precisa sobre las variaciones espacio-temporales que ocurren en los ecosistemas a partir de las bandas espectrales de una imagen, dado que asumen una relación (empírica o teórica) con variables bio-geofísicas de la superficie.

El caso más conocido es el Índice de Vegetación de Diferencia Normalizada (NDVI, Normalized Difference Vegetation Index) (Tucker, C. J. 1979). Las bases teóricas para su desarrollo se derivan del comportamiento de la firma espectral típica de la vegetación. La energía reflejada en el visible es muy baja debido a la actividad de los pigmentos fotosintéticos que determinan una fuerte absorción en las porciones del espectro correspondientes al azul (470 nm) y al rojo (670 nm). Por el contrario, casi toda la radiación en el infrarrojo cercano se refleja (hay muy poca absorción) dependiendo del índice de área foliar (LAI), la distribución angular de las hojas y su anatomía y morfología. El contraste entre las respuestas del rojo (R) y el infrarrojo cercano (IRc) constituye una medida sensible de la cantidad de vegetación, donde la máxima diferencia entre R y IRc corresponde al estadío de mayor densidad o vigor de la vegetación, y el mínimo contraste a áreas de muy poca vegetación o vegetación en estado senescente. Este salto se denomina “red edge” o “borde rojo”

El más sencillo de estos índices es simplemente el cociente entre las bandas roja e infrarroja:

$ VI= IRc/ R$

El NDVI, mencionado anteriormente corresponde al cociente:

$ NDVI = (nir -red) / nir +red)$

Como cociente normalizado tiene ciertas ventajas. Reduce parcialmente los efectos ambientales (atmósfera) que se observan en las bandas individuales, variaciones provocadas por la topografía, sombras, variaciones en iluminación. Realza variaciones pequeñas en las respuestas espectrales de las coberturas. Finalmente, este índice se puede determinar a partir del contaje digital crudo, radiancias, reflectancias TOA, reflectancias de superficie. Sin embargo, si bien las unidades se cancelan, y el cociente neutraliza relativamente el ruido introducido por la atmósfera, los resultados indican que es recomendable trabajar con las imágenes en reflectancia de superficie (corregidas atmosféricamente), especialmente cuando se trabaje con imágenes de más de una fecha (o con series temporales).

En algunos casos puede ser interesante reemplazar la banda del rojo visible por la del verde visible, donde entrarían en juego otros pigmentos y la cual es sensible al estrés en las plantas.

El NDVI describe adecuadamente el comportamiento de la vegetación en ambientes donde la cobertura vegetal es alta o por lo menos hay baja proporción de suelo desnudo. En casos donde la cobertura de la vegetación es muy pobre una buena parte de la reflectancia registrada por los sensores corresponde a la contribución del suelo y el NDVI no tiene un comportamiento eficiente frente a estas situaciones (Price, J. C. 1987). El Índice verde corregido por suelo (SAVI, Soil Adjusted Vegetation Index), da cuenta de este fenómeno (Huete, A. R. 1988):

$ SAVI= ((nir -red)/(nir +red+L))*(1+L)$

Donde L es la contribución del suelo y puede variar entre 0 y 1, según el porcentaje de suelo desnudo. Un valor de 0,5 indicaría una proporción de 50% de suelo desnudo.

El objetivo en este caso es analizar la potencial información que aporta el uso de diferentes índices espectrales y evaluar la información que recuperan sobre las características biofísicas de los ambientes y su variación temporal.

Se utilizará la imagen con reflectancia en superficie de Landsat y vamos a cálcular varios indices. También puede usar otra más que le parezca interesante.

NDVI

vi <- function(img, k, i) {
  bk <- img[[k]]
  bi <- img[[i]]
  vi <- (bk - bi) / (bk + bi)
  return(vi)
}


ndvi <- vi(imagen_stack, 5, 4)
plot(ndvi, col=rev(terrain.colors(10)), main = "Chaco-NDVI in 2022")

EVI

vi2 <- function(img, m, k, i) {
  bm <- img[[m]]
  bk <- img[[k]]
  bi <- img[[i]]
  vi2 <- 2.5 * ((bm - bk) / (bm + 6* bk - 7.5* bi + 1))
  return(vi2)
}

# For Landsat 8: NIR = 5, red = 4, blue=2.
evi <- vi2(imagen_stack, 5, 4, 2)
plot(evi, col=rev(terrain.colors(6)), main = "Montes- EVI in 2013")

par(mfrow=c(1,2))
hist(ndvi, main="NDVI")
hist(evi, main="EVI")

# correlation plot between indices

indices_ndvi<- as.data.frame(ndvi, xy=TRUE, na.rm=TRUE) # convertir a data frame
indices_evi<- as.data.frame(evi, xy=TRUE, na.rm=TRUE) # convertir a data frame
indices <- merge(indices_ndvi,indices_evi,by=c("x","y"))
names(indices) <- c('x','y',"ndvi", "evi")
ggplot(indices, aes(x=evi, y=ndvi)) + 
  geom_point()+
  geom_smooth(method=lm)
## `geom_smooth()` using formula 'y ~ x'

#plot(indices_ndvi$layer,indices_evi$layer)

Índice Diferencial de Agua Normalizado (NDWI)

El NDWI se utiliza como una medida de la cantidad de agua que posee la vegetación o el nivel de saturación de humedad que posee el suelo. Para obtener este índice la combinación de bandas es la siguiente: Landsat 9 (3-6)/(3+6)

wi <- function(img, j, k) {
  bj <- img[[j]]
  bk <- img[[k]]
  wi <-  (bj - bk) / (bj +bk)
  return(wi)
}
ndwi <- wi(imagen_stack, 3,5)
plot(ndwi, col=rev(terrain.colors(10)), main = "NDWI")